{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "sealed-vertical",
   "metadata": {},
   "source": [
    "The purpose of this notebook is to evaluate how many query vectors (nq) we need \n",
    "to get an accurate estimate of the intersection @ 10 measure. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "id": "hairy-relaxation",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import faiss\n",
    "from faiss.contrib import datasets\n",
    "from matplotlib import pyplot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "id": "emotional-filter",
   "metadata": {},
   "outputs": [],
   "source": [
    "# start with BigANN 10M\n",
    "\n",
    "ds = datasets.DatasetBigANN(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "alpha-webcam",
   "metadata": {},
   "outputs": [],
   "source": [
    "# we need more query vectors than the usual 10k, so pick them from the training set\n",
    "\n",
    "xtt = ds.get_train(maxtrain=2 * 10**6)\n",
    "\n",
    "big_xq = xtt[:10**6]  # 1M queries \n",
    "xt = xtt[10**6:]      # 1M training vectors"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "imported-variance",
   "metadata": {},
   "source": [
    "# Ground truth"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "adequate-jefferson",
   "metadata": {},
   "outputs": [],
   "source": [
    "index = faiss.IndexFlatL2(128)\n",
    "index.add(ds.get_database())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "minus-defendant",
   "metadata": {},
   "outputs": [],
   "source": [
    "index = faiss.index_cpu_to_all_gpus(index, ngpu=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "trying-andorra",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 4min 59s, sys: 2min 8s, total: 7min 8s\n",
      "Wall time: 7min 8s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "Dgt, Igt = index.search(big_xq, 10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "maritime-oriental",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.save(\"/tmp/Dgt.npy\", Dgt)\n",
    "np.save(\"/tmp/Igt.npy\", Igt)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "supreme-jewel",
   "metadata": {},
   "source": [
    "# With some index"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "comic-therapy",
   "metadata": {},
   "outputs": [],
   "source": [
    "# we pick some index that has one search-time parameter \n",
    "# to represent different speed-accuracy tradeoffs.\n",
    "\n",
    "index = faiss.index_factory(128, \"IVF16384,SQ4\")\n",
    "index = faiss.index_cpu_to_all_gpus(index, ngpu=1)\n",
    "\n",
    "index.train(xt)\n",
    "\n",
    "index.add(ds.get_database())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "id": "perfect-associate",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1\n",
      "4\n",
      "16\n",
      "64\n"
     ]
    }
   ],
   "source": [
    "# run queries \n",
    "\n",
    "res_per_nprobe = {}\n",
    "for nprobe in 1, 4, 16, 64:\n",
    "    print(nprobe)\n",
    "    D, I = [], []\n",
    "    index.nprobe = nprobe\n",
    "    # Faiss crashes when searching all at once\n",
    "    for i0 in range(0, 10**6, 10000): \n",
    "        Di, Ii = index.search(big_xq[i0 : i0 + 10000], 10)\n",
    "        D.append(Di)\n",
    "        I.append(Ii)\n",
    "    D = np.vstack(D)\n",
    "    I = np.vstack(I)\n",
    "    res_per_nprobe[nprobe] = I"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "novel-yukon",
   "metadata": {},
   "source": [
    "# Stats on intersection measure"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "id": "polish-sheep",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "nprobe     1 inter 25.914 %\n",
      "nprobe     4 inter 49.756 %\n",
      "nprobe    16 inter 68.647 %\n",
      "nprobe    64 inter 77.002 %\n"
     ]
    }
   ],
   "source": [
    "# evaluate intersection measures for the 1M queries\n",
    "\n",
    "for nprobe in 1, 4, 16, 64: \n",
    "    I = res_per_nprobe[nprobe]\n",
    "    ninter = faiss.eval_intersection(I, Igt)\n",
    "    print(f\"nprobe {nprobe:-5d} inter {100 * ninter/I.size:.3f} %\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "id": "consistent-service",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZgAAAEGCAYAAABYV4NmAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAB1OUlEQVR4nO2dZ3hVxdaA35XeEwghhASk14B0FUIRlCIioIgdEBS7YL3y3atiL6gU7xWlIyqIIAjSFYL03pt0CQQCCamkZ74feyckIT3n5CRh3ufZz9ln7T0za3JOztozs2YtUUqh0Wg0Go2lsbO1AhqNRqOpnGgDo9FoNBqroA2MRqPRaKyCNjAajUajsQrawGg0Go3GKjjYWoHyQrVq1VSdOnVKXD4hIQF3d3fLKVQBuNn6fLP1F3SfbxZK0+ddu3ZdUUr55XVNGxiTOnXqsHPnzhKXDw0NpVu3bpZTqAJws/X5Zusv6D7fLJSmzyJyNr9reopMo9FoNFbBagZGRFxEZLuI7BORQyLynikfKyLnRWSvedyTrcwYETkhIsdEpFc2eVsROWBemyQiYsqdReRnU75NROpkKzNURI6bx1Br9VOj0Wg0eWPNKbJkoLtSKl5EHIGNIrLCvDZeKfVF9ptFpBnwMNAcqAn8ISKNlFLpwGRgJLAVWA70BlYAI4CrSqkGIvIw8BnwkIhUBd4F2gEK2CUiS5RSV63YX41Go9Fkw2oGRhkxaOLNt47mUVBcmv7APKVUMnBaRE4AHUTkDOCllNoCICLfAwMwDEx/YKxZfgHwX3N00wtYo5SKMsuswTBKc4vTh9TUVMLCwkhKSir0Xm9vb44cOVKc6is8Je2zi4sLQUFBODo6WkErjUZTXrDqIr+I2AO7gAbA/5RS20SkD/CiiAwBdgKvmSOLQIwRSiZhpizVPM8tx3w9B6CUShORGMA3uzyPMkUmLCwMT09P6tSpgzkrly9xcXF4enoWt4kKTUn6rJQiMjKSsLAw6tatayXNNBpNecCqBsac3molIj7AIhEJxpju+gBjNPMB8CUwHMjrF1wVIKeEZbIQkZEYU2/4+/sTGhqa47q3tze+vr7Ex8fnLnoD6enpxMXFFXpfZaKkfXZyciI6OvqGv3d5Jz4+vsLpXFpupj67bNuOx2+/UT0qioNVqxLfvz9Jt3WwtVplgrU+5zJxU1ZKRYtIKNA7+9qLiEwFfjffhgG1shULAi6Y8qA85NnLhImIA+ANRJnybrnKhOah1xRgCkC7du1Ubje9I0eO4OXlVaQ+6hFM8XBxcaF169YW1si6aPfVykvM0qWEz52LMqfD7aOi8Jk7l4BmTfHu18/G2lkfa33O1vQi8zNHLoiIK3AXcFREArLdNhA4aJ4vAR42PcPqAg2B7UqpcCBORG4311eGAL9lK5PpITYIWGuu/awCeopIFRGpAvQ0ZRqNRnMDEeMnZBmXTFRSEhHjJ9hGoUqCNffBBADrRGQ/sANj0f134HPT5Xg/cCfwCoBS6hAwHzgMrAReMKfYAJ4DpgEngJMYC/wA0wFf0yHgVeAts64ojOm3HebxfuaCv8Zg2LBhLFiwoMTlhw8fTr169QgODragVhqNbUgLDy+WXFM0rOlFth+4YQ5EKfVEAWU+Aj7KQ74TuOGXTCmVBDyYT10zgBnFULnULN5znnGrjnEhOpGaPq680asxA1oX27fAYqSlpeHgYJ2PeNiwYTz55JM899xzVqlfoylLHAICSLtwIU+5puTonfwWYtnBS4z59QDnoxNRwPnoRMb8eoDFe86XuM4zZ87QtGlTnn76aZo3b07Pnj1JTEwEoFu3bowePZqOHTsSHBzM9u3bARg7diwjR46kZ8+eDBkyhLNnz9KjRw9atmxJjx49+Oeff7Lq/+OPP+jcuTONGjXi99+NpbD09HTeeOMN2rdvT8uWLfnuu+/y1K1Lly5UqVKlxH3TaMoTfi+/DLk8RcXFheqvjLaNQpUEHYusiLy39BCHL8Tme33PP1dJSc/pqJaYms6bC/Yzd/s/eZZpVtOLd/s1L7Dd48ePM3fuXKZOncrgwYNZuHAhjz/+OGAEqNu8eTN//fUXw4cP5+BBYzlr165dbNy4EVdXV/r168eQIUMYOnQoM2bM4OWXX2bx4sWAYcDWr1/PyZMnufPOOzlx4gTff/893t7e7Nixg+TkZDp16kTPnj21S7GmcpOWCkphX6UKaVevIoBH9+43xQK/NdEjGAuR27hcl2eUqt66devSqlUrANq2bcuZM2eyrj3yyCOAMZqIjY0lOjoagPvuuw9XV1cAtmzZwqOPPgrAE088wcaNG7PKDx48GDs7Oxo2bEi9evU4evQoq1ev5vvvv6dVq1bcdtttREZGcvz48VL1QaMpz2SkpHD5m29wadmShps3EfHtZNw73kHizp2olBRbq1eh0SOYIlLYSOOOj/8gPDb5Bnmgjys/P3NHidt1dnbOOre3t8+aIgNu2PyZ+b6gsNvZy+RVXinF119/Ta9evXIX1WgqJdE/zyftQjgBH3yQ9T9RdfgIzj31FDG/L8Pn/oE21rDiokcwFmLUnXVwdbTPIXN1tOeNXo2t1ubPP/8MwMaNG/H29sbb2/uGezp27Mi8efMA+PHHHwkJCcm69ssvv5CRkcHJkyc5deoUjRs3plevXkyePJnU1FQA/v77bxISEqzWB43GlmRcu8aV777DrUMH3Dt2zJK7d+qIc+PGRM2cgbHzQVMS9AjGQvQN9sfFxbVMvciqVKlCx44diY2NZcaMvB3mJk2axPDhwxk3bhx+fn7MnDkz61rjxo3p2rUrly5d4ttvv8XFxYWnnnqKM2fO0KZNG5RS+Pn5Za3ZZOeRRx5h3bp1REZGEhQUxHvvvceIESOs1VWNxipE/fgj6Veu4Ddp0g2je9/hT3LhX2+RsGEDHl262FDLCoxSSh9K0bZtW5Wbw4cP3yDLj9jY2CLfawm6du2qduzYUaZt5qY0fS7O37a8sG7dOlurUOZU5j6nxcaqox1uU2dHjswhz+xzRnKy+rtLV3Vm6DAbaFe2lOZzBnaqfH5X9RSZRqO5KYmaOZOMmBiqjxqV53VxcqLqkCFc27qVxEOHyli7yoE2MBWU0NBQ2rVrZ2s1NJoKSVpUFFGzZuPZqxcuzZrle5/P4Aexc3cnasbMfO/R5I82MBqN5qYjcuo0MpKS8Hv5pQLvs/f0xOehh4hduZLU8yXfNH2zog2MRqO5qUi9FMHVn37C+777cK5fv9D7qz7xOIgQ9f33ZaBd5UIbGI1Gc1Nx5dvJqIwMqr34QpHudwwIwLvvPVz9ZQHpMTFW1q5yoQ2MRqO5aUg5d47oXxbgM+gBnIKCCi9gUnX4cNS1a1yd97MVtat8aANzk1LacP1gBMZs3bo19957r4W00misy5X//g+xt6fas8WLAu7SuDHunToR9cMcMnT4mCKjDYwl2T8fxgfDWB/jdf98m6qTlpZm1fonT55M06ZNrdqGRmMpkk+eJGbpUqo8+iiO/tWLXd53xHDSL18hdulSK2hXOdEGxkI4HFkES1+GmHOAMl6XvlwqI1Oew/WHhYWxatUqnnrqqRL3T6MpSy5P+ho7Fxd8Rz5dovJud9yBc9OmRM6YicooXRDbmwUdKqaorHgLLh7I97JL2HZIzzV0Tk2E316EXbPzLlSjBfT5tMBmy2u4/tGjR/P++++Tof/RNBWApMOHiVu1imrPP49DCfMYZYWPeeNN4v/6C08r5LCvbOgRjKXIbVyy5DdGWC4O5TFc/++//0716tVp3fqGhKUaTbkkYuJE7Ly9qfrksFLV49W7Nw4BAXrjZRHRI5iiUshIQ33ZDInLYyOWdy14clmJmy2P4fo3bdrEkiVLWLZsGcnJycTGxvL444/zww8/FK1TGk0Zcm33bhLW/4Xfa69i7+lZqrrE0ZGqQ4YQ8dlnJB44iGuLGzK5a7KhRzAWIrnzW+DomlPo6Ao93rFam7YK1//JJ58QFhbGwYMHmTdvHt27d9fGRVMuUUpxefwE7KtVo+pjj1mkTp8HB2Hn4UHUzLwjmGuuYzUDIyIuIrJdRPaJyCERec+UVxWRNSJy3Hytkq3MGBE5ISLHRKRXNnlbETlgXpsk5qO3iDiLyM+mfJuI1MlWZqjZxnERGWqtfmaS1nQg9JtkjFgQ47XfJGg52GptZobrf/bZZ5k+fXqe90yaNImZM2fSsmVL5syZw8SJE7OuZYbr79OnT45w/c2aNaNNmzYEBwfzzDPPWN0bTaOxFgmbN3Ntxw6qPfssdm5uFqnT3sODKg8/ROzKVaSEhVmkzkpLfmGWS3uAkdbaPHcEtgG3A58Db5nyt4DPzPNmwD7AGagLnATszWvbgTvMOlcAfUz588C35vnDwM/meVXglPlaxTyvUpC+Olx/8dHh+is/FbnPGRkZ6tSgB9Xfd96p0pOTi1yuKH1OuXhRHQ5uocI/+LAUGpYfKly4frPtePOto3kooD+Q6VY1GxhgnvcH5imlkpVSp4ETQAcRCQC8lFJbzM58n6tMZl0LgB7m6KYXsEYpFaWUugqsAXpbp6cajaY8Er92LUkHDuD3wgvYOTlZtG5Hf3+8+/YleuFC0k3nGs2NWHWRX0TsgV1AA+B/SqltIuKvlAoHUEqFi0jmjqdAYGu24mGmLNU8zy3PLHPOrCtNRGIA3+zyPMpk128kMBLA39+f0NDQHNe9vb2Ji4srUl/T09OLfK8lWGpu9irLNnNTmj4nJSXd8Pcu78THx1c4nUtLhe1zRgZVP/oY8a/OXh8fKEYfitpn+xbBVFu8mN0ff0LCPX1KrGp5wFqfs1UNjFIqHWglIj7AIhEpyOVC8pCpAuQlLZNdvynAFIB27dqpbrn82o8cOYJnEb1O4uLiinxvZaE0fXZxcalwbs6hoaHk/o5Udipqn2OW/s6FCxeo+eUXtOjRo1hli9Pnf0LXY795M20/eB+7bB6fFQ1rfc5l4kWmlIoGQjGmqS6Z016YrxHmbWFArWzFgoALpjwoD3mOMiLiAHgDUQXUpdFoKjkqNZXL//0a58aN8epj3ZGF74jhpF+5QsySJVZtp6JiTS8yP3Pkgoi4AncBR4ElQKZX11DgN/N8CfCw6RlWF2gIbDen0+JE5HZzfWVIrjKZdQ0C1prrNKuAniJSxfRS62nKNBpNJSd68WJSz/6D36hRiJ11n6HdbrsN52ZNiZo5S4ePyQNr/vUDgHUish/YgbHo/jvwKXC3iBwH7jbfo5Q6BMwHDgMrgRfMKTaA54BpGAv/JzE8yQCmA74icgJ4FcMrDaVUFPCB2e4O4H1TptFoKjEZKSlc+WYyLre2xOPOblZvzwgfM4KUU6eID11v9fYqGtb0ItuvlGqtlGqplApWSr1vyiOVUj2UUg3N16hsZT5SStVXSjVWSq3IJt9p1lFfKfWiOUpBKZWklHpQKdVAKdVBKXUqW5kZpryBUkrHdchFacP1Dx8+nHr16hEcfOOy2tdff03jxo1p3rw5b775ZmnU1GiKRfS8n0kLD6f66NE3RKqwFl69euJQM4CoGXrjZW70Tn4LsuzUMnou6EnL2S3puaAny06VPESMJbDmBslhw4bx66+/3iBft24dv/32G/v37+fQoUO8/vrrVtNBo8lOxrVrXPnuO9xuuw33O+4os3bF0RHfoUO5tnMnifv3l1m7FQFtYCzEqn9WMXbzWMITwlEowhPCGbt5bKmMTHkO19+lSxeq5BGVdvLkybz11ltZMdSqVy9+3g2NpiREzfmB9MhI/EaPKvO2vR8YhJ2nJ5HT9SgmOzrYZRH5bPtnHI06mu/1fZf3kZqRmkOWlJ7EO5veYcHfeU9FNanahH91+FeB7ZbXcP358ffff7Nhwwb+/e9/4+LiwhdffEH79u2LVFajKSnpsbFETp+OR7duuNnA/d3ew50qDz9M5PTppPzzD061a5e5DuURPYKxELmNSyYpGaVLr1oew/UXRFpaGlevXmXr1q2MGzeOwYMHZ4YO0misRuTMmWTExuI36mWb6VDl8cfB3p6oWfnkf7oJ0SOYIlLYSOOu+XdxKfHSDfIA9wBm9i65j0F5DNdfEEFBQdx///2ICB06dMDOzo4rV67g5+dXovo0msJIi4ri6uzv8ezTGxcbpvB29K+Od79+RP/6K9VeerHEic0qE3oEYyGebf4sLvYuOWQu9i6MamO9+WBbhesviAEDBrB27dqssikpKVSrVq3EfdRoCiNyylQykpLwe+klW6uC75PDUElJXJ0719aqlAv0CMZC9KrdC1dXVybunsjFhIvUcK/BqDaj6Fuvr9XazAzXHxsby4x8XCQnTZrE8OHDGTduHH5+fsyceX00lRmu/9KlSznC9Z85c4Y2bdqglMLPzy9rzSY7jzzyCOvWrSMyMpKgoCDee+89RowYwfDhwxk+fDjBwcE4OTkxe/bsMnMX1dx8pF68yNWffsK7f3+c69WztTo4N2yIe9cuXP3hR3xHjKjQ4WMsQn5hlm+2Q4frLz46XH/lp7z3+cI776rDwS1U8rkwi9VZ2j7Hb92mDjduoqLm/WwZhcqACheuX6PRaKxJyrlzRC9cSJUHH8Qp6IZg6TbDrUN7XJo3J2rmzJs+fIw2MBWU0NBQ2rVrZ2s1NBqbceW//0UcHPB99hlbq5IDEcF3xHBSzpwhft06W6tjU7SB0Wg0FY7k48eJWbKUKo89imM53Mzr2bMnjoGBN/3GS21gNBpNhePypK+xc3PD96mnbK1KnoiDA1WHDiVx926u7dlja3VshjYwGo2mQpF48BBxa9ZQddiwcr3XxOeB+7Hz9iZqxs0ba1cbGI1GU6G4PHEi9t7eVH1ymK1VKRA7dyN8TNwff5By9qyt1bEJ2sDcpJQ2XH90dDRPPPEETZo0oWnTpmzZsiXH9S+++AIR4cqVK6VVVaPJ4tquXSRs2IDvyKex9/CwtTqFUuWxRxEHByJnzbK1KjZBGxgLErN0Kce79+BI02Yc796DmKVLbaqPNcP1jxo1irvuuoujR4+yb98+mmYL0XHu3DnWrFlDbR3wT2NBlFJcHj8Be79qVDHj65V3HKtXx6v/fcT8uoi0qJsv56E2MBYiYcVKwt9+h7QLF0Ap0i5cIPztd0plZMpruP7Y2Fj++usvhgwZAoCTkxM+Pj5Z11955RU+//xzvYNfY1ESNm3m2s6dVHv2WezMYK4VAd8nn0QlJ3P1p5svfIwOFVNELn78MclH8g/Xf23vXkjNGVFZJSUR/u//ED3/lzzLODdtQo3/+78C2y2P4fpPnTqFn58fzz33HIcPH6Zt27ZMnDgRd3d3lixZQmBgILfeemthf1KNpsgopbg8YQKONWtS5cEHba1OsXCuXx+Pbt24+uOP+I4YXqGMY2nRIxhLkZp3uH6VUvnC9aelpbF7925GjBjBnj17cHd359NPP+XatWt89NFHvP/++6Xqs0aTm7g//iDp4EGqvfgi4uRka3WKje+I4aRfvUrMb7/ZWpUyRY9gikhhI42/u91J+sWLN8gdatbkljnfl7jd8hiuPygoiKCgoKxEYoMGDeLTTz/l5MmTnD59Omv0EhYWRps2bdi+fTs1atQorKsaTZ6o9HSuTJqEU926eN/Xz9bqlAjXdu1wadGCyJkz8XnwQcTe3tYqlQlWG8GISC0RWSciR0TkkIiMMuVjReS8iOw1j3uylRkjIidE5JiI9MombysiB8xrk8T8ZRQRZxH52ZRvE5E62coMFZHj5jHUWv3MxOv55xGXnOH6xcWF6q+MtlqbtgrXX6NGDWrVqpU1svnzzz9p1qwZLVq0ICIigjNnznDmzBmCgoLYvXu3Ni6aUhG7fDnJx0/g9/JLiEPFfCbODB+TevYf4sx0FjcD1vy00oDXlFK7RcQT2CUia8xr45VSX2S/WUSaAQ8DzYGawB8i0kgplQ5MBkYCW4HlQG9gBTACuKqUaiAiDwOfAQ+JSFXgXaAdoMy2lyilrlqrs+59euPq6kLE+AmkhYfjEBBA9VdG493Pek9ctgzX//XXXzN8+HDS0tKoV69ejno1GkuhUlO5/PV/cW7SBM8SJsErL3jedReOQUFETZ+B191321qdsiG/MMuWPoDfgLuBscDreVwfA4zJ9n4VcAcQABzNJn8E+C77Pea5A3AFkOz3mNe+Ax4pSD8drr/46HD9lR9b9zlq3s/qcOMmKnbt2jJr05p9jpzzgzrcuIlK2LXbam2UBGuF6y+T8aY5ddUa2AZ0Al4UkSHAToxRzlUgEGOEkkmYKUs1z3PLMV/PASil0kQkBvDNLs+jTHa9RmKMjPD39yc0NDTHdW9vb+Li4orUx/T09CLfawnS09NJSEgo0zbz0qGk7SclJd3w9y7vxMfHVzidS4tN+5yaSrUJE0ivW5ddAGWkh1X77FcNP3d3jn32GTHPPWudNkqAtfpsdQMjIh7AQmC0UipWRCYDH2BMXX0AfAkMxxh55EYVIKeEZa4LlJoCTAFo166d6tatW47rR44cwdPTM4+qbiQuLq7I91qCDRs2lFlb+VGaPru4uNC6dWsLa2RdQkNDyf0dqezYss9Rs2dz6epV6o4fj/vtt1m9vWWnljFx90TCE8IJSAqwWkbaiKNHifz2O5rdcgvO2dz/bYm1PmeruimLiCOGcflRKfUrgFLqklIqXSmVAUwFOpi3hwG1shUPAi6Y8qA85DnKiIgD4A1EFVBXsTFGgBpLov+mmoKIWbqU493u5NInnyJOTqRdjrB6m8tOLWPs5rGEJ4QDEJ4QztjNY1l2apnF26r62GOIoyNRs2ZbvO7yhtVGMKan13TgiFLqq2zyAKVUuPl2IHDQPF8C/CQiX2Es8jcEtiul0kUkTkRux5hiGwJ8na3MUGALMAhYq5RSIrIK+FhEMkOt9sRY4ykWLi4uREZG4uvrq3elWwilFJGRkbjk8rjTaMAwLmH/+Td2yYYXo0pJIew//wYoscNMekY68anxxCbHEpsSS0xKDLEpscSlxGXJ5h2dR1J6Uo5ySelJTNw90eKjGIdq1fDu35+YRYvwe/klHHx9LVp/ecKaU2SdgCeAAyKy15T9H/CIiLTCmLI6AzwDoJQ6JCLzgcMYHmgvKMODDOA5YBbgiuE9tsKUTwfmiMgJjJHLw2ZdUSLyAbDDvO99pVSxAwEFBQURFhbG5cuXC703KSnppvvRLGmfXVxcCAoKKvxGzU3H2XEf45icc9OyXXIqZz//mKC7O2UZhNiU2KzzuJS4G2SZ53EpccSlFrxO6GjnSGpG3hulLybcuLfNElR98kmif/mFqz/+hN/LL1mljfKA1QyMUmojea+FLC+gzEfAR3nIdwLBeciTgDzjRiilZgClSifn6OiYI0RKQYSGhla4NYXScjP2WWNdHCKi85TbX46m689d8y3nYu+Cl5MXXs5eeDl5UcOtBg19Gma9z34t8/B08sTL2QsXexd6LeyVNT2WnRru1tnD5VyvLh7du3P1p5/wffqpShs+pmLuWtJoNJWSK17gF3ujPNILxnQYk6/BcLIvXfiYUW1GMXbz2BzTZE72ToxqM6pU9RaE74jhnF27luhFi6haQaJDFxdtYDQaTbkh9DZXHlyTmEOW5AArelblk6bW+xHOXGfJ9CIThAbeDaziRZaJa5s2uNzakqhZs6ny0EOVMnyMDnap0WjKBcevHsc1NoV04IonZACXvWDmvc6EDC84FqAl6FuvL6sHrebrW75mZMuRHI46zImrJ6zWnojgO3wEqf/8Q9wff1qtHVuiDYxGo7E5CakJvP7nK3Q5qIi7vQkf/KsWj4xx5MM3a9Fr5EdWHUnkxeNNH8fVwZVpB6dZtR3Pu3rgWLs2kTOmV0r3fT1FptFobIpSivc2v4ffnrN4JWQQNGwUq228odXHxYfBjQYz58gcXrj1BWp51Sq8UAkQe3uqDhvKpfc/4HjnzqRHRpVJHMOyQo9gNBqNTfnl719YcWYFT56qhUP16nhki/htS4Y2H4qDODDjUKmcUQtFzJQc6VciLZYNt7ygDYxGo7EZhyIP8en2T+np3o6q+87iPXBguQnJ7+fmx8CGA/ntxG9cSrhktXau/O+bG2QqKYmI8ROs1mZZoQ2MRqOxCbEpsbwW+hpVXaoyOrwlZGTg88D9tlYrB8OaDyNDZTDr0CyrtZEWfuP+m4LkFQltYDQaTZmjlOLtjW9zKeESX3T+nOQlK3C77Tacate2tWo5CPIMom+9viz4ewFRScUOBlIkHAICiiWvSGgDo9FoypzvD3/P2nNreaXtKzQ8k0LquXP4DHrA1mrlyYgWI0hOT+aHwz9Ypf7qr4wu82y4ZYU2MBqNpkzZG7GXCbsm0KN2D55o9gTRCxZi5+WFZznN8ljPux533XIXc4/OJTYljzADpcS7Xz8CPngfezPopX2VKgR88L72ItPYjsV7ztPp07XUfWsZnT5dy+I9522tkkZTKFeTrvL6+tep4V6D9zu9T0ZsLHGrV+N9773YleNgsU+3eJr41Hh+PvqzVer37tePhqHrsHN3x/PuuyuFcQFtYCoki/ecZ8yvBzgfnYgCzkcnMubXA9rIaMo1GSqDMRvHEJUUxZfdvsTLyYuYpb+jUlLK7fRYJk19m9I5sDNzDs/hWuo1q7Qhjo643X47CRs3VppNl9rAVEDGrTpGYmp6DlliajrjVh2zkUYaTeFMOzCNTec38VaHt2jm2wyA6IULcW7WFJdmzWysXeE83fJpriZfZeHxhVZrwyOkE6kXLpBy+ozV2ihLCjUwIlJfRJzN824i8rKI+FhdM02eKKU4H52Y57UL+cg1GluzPXw7/9v7P/rU7cODjYwMG4mHDpF85Ag+gwbZWLui0bp6a9r5t2PWwVmkpKdYpQ13c5NpwsaNVqm/rCnKCGYhkC4iDTASfNUFfrKqVpo8uRiTxJAZ2/O9bm8nnIgoOLmSRlPWXEm8wpt/vcktXrcw9o6xWdlhYxYuRJyd8b73XhtrWHSebvk0EYkR/HbyN6vU71SrFk633EL8ppvHwGQopdIw0htPUEq9AlR8B+0KxtJ9F+g14S92nrnKoLZBuDrm/Oic7O1wdrCj76SNfL/lTKWZw9VUbNIy0njzrzdJSE3gy65f4uboBkBGUhIxS3/Hs2dP7L28bKxl0bkj4A6CfYOZcWAGaRlpVmnDPSSEa9t3kJFinVFSWVIUA5MqIo8AQ4HfTZmj9VTSZCfmWiovz93DS3P3ULeaO8tHdeaLB2/lk/tbEujjigCBPq58Pqgl697oxh31fXnnt0M8OWsHEXFJhdav0ViTb/Z+w46LO/jP7f+hYZWGWfK41avJiIurMNNjmYgIT7d8mrD4MFaeWWmVNtxDOqESE0nctcsq9ZclRQn68yTwLPCRUuq0iNQFrLPjSJODDccv88Yv+7kSn8xrdzfiuW71cbA3ngkGtA5kQOvAG8rMHNaeH7ae5cNlR+g9YQOf3t+Cns2tk/ZVoymIDWEbmHpgKvc3vJ/+DfrnuBa9YCGOtWvj1qG9jbQrOd1qdaOBTwOm7Z/GPXXvwU4s6yvl3qEDODoSv2Ej7nfcYdG6y5pC/zJKqcPAv4Dd5vvTSqlPra3YzUxiSjrv/naQJ6Zvx8PFgUXPd+KlHg2zjEtBiAhP3FGHZS+HEODtwsg5uxjz636upVhnOK/R5MXFhIv838b/o1GVRozpMCbHtZSzZ7m2fTs+99+ftR5TkbATO55q8RQnY06y7p91lq/f3R23Nm0qxUJ/UbzI+gF7gZXm+1YisqQI5WqJyDoROSIih0RklCmvKiJrROS4+VolW5kxInJCRI6JSK9s8rYicsC8NknMb6WIOIvIz6Z8m4jUyVZmqNnGcREZWvQ/iW3Zdy6avl9vYPaWswzvVJffXwqhRZB3setpUN2TRc934tmu9Zm34xx9J21k77loyyus0eQiNT2V19a/Rkp6Cl92/RIXh5wbKKMX/gp2dngPHGAbBS1Arzq9qOVZiykHplhlvdM9pBPJf/9N6qUIi9ddlhRlbDcW6ABEAyil9mJ4khVGGvCaUqopcDvwgog0A94C/lRKNQT+NN9jXnsYaA70Br4Rkcwk1ZOBkUBD8+htykcAV5VSDYDxwGdmXVWBd4HbTN3fzW7IyiOp6RmMX/M390/eTFJKOj89dRvv9GuGi2PJ83Q7OdjxVp8mzH36dpJT03lg8ma+/vM4aekZFtRco8nJhN0T2H95P+91eo863nVyXFNpacQsWoRHly44+vvbRkEL4GDnwIjgERyOPMzmC5stXr9H584AJGzaZPG6y5KiGJg0pVRMLlmhJlspFa6UypxWiwOOAIFAf2C2edtsYIB53h+Yp5RKVkqdBk4AHUQkAPBSSm1RxqPC97nKZNa1AOhhjm56AWuUUlFKqavAGq4bpXLHycvxDJq8mYl/Hue+W2uyYnQXOjaoZrH6b6/ny4rRXejbIoAv1/zNw1O2ci7KOruRNTc3f579k+8Pf88jTR6hd50b/+XiN2wg7fLlcr9zvyjcV/8+/N38mbJ/isXrdm7cGHu/ahV+mqwoi/wHReRRwF5EGgIvA8Uy2ebUVWtgG+CvlAoHwwiJSHXztkBga7ZiYaYs1TzPLc8sc86sK01EYgDf7PI8ymTXayTGyAh/f39CQ0OL060cxMfHF7t8hlKs/SeN+cdScLSHF1o5094/mj3brPPUcn8ABChnvj98lbu/XMfjzZzoVNOhxPPgJelzReZm6y8Ur89XUq/wefjn1HaqTftr7fMs5/3tdzh6ebFLBMrp37I4fQ5xDmFhxEKmrZxGA5cGFtXDq359Utev5/jatWBn3aAr1vpuF8XAvAT8G0gG5gKrgA+K2oCIeGBs1hytlIot4McsrwuqAHlJy1wXKDUFmALQrl071a0UecBDQ0MpTvmLMUm8sWAfG45foVtjPz5/oCXVvawf7K8b8MTVa7w6fx/TDkQRji8fDQzGx82p2HUVt88VnZutv1D0PienJ/PE8idwcHTgu3u/I8gz6IZ70i5f5vjBg/g+OYzgHj2soK1lKM7nfFvabaxbuI6dDjt5qttTFtUjJi6eC1u3cZufH64tWli07txY67tdFC+ya0qpfyul2iul2pnnRdpgISKOGMblR6XUr6b4kjnthfmauYoVBtTKVjwIuGDKg/KQ5ygjIg6ANxBVQF3lgt/2nqfn+PXsPHOVjwYGM3NY+zIxLpkEVXFj7tO382bvxqw6dJHeEzaw6cSVMmtfU/kYt2McR6KO8FGnj/I0LgAxv/0G6el431/xp8cycXVw5YlmT7Dp/CYORR6yaN3unTqCCPEbNli03rIkXwMjIhPM16UisiT3UVjF5lrIdOCIUuqrbJeWYGzaxHz9LZv8YdMzrC7GYv52czotTkRuN+sckqtMZl2DgLXmOs0qoKeIVDEX93uaMpsSfS2FF3/azah5e2lQ3YMVozrz2G232MRV095OeL5bAxa/0Ak3Z3sem7aND38/TFKuIJoaTWEsO7WMn4/9zJPNn+TO2nfmeY9SiugFC3Ft2xbnekXxEao4PNT4ITwdPZm2f5pF63WoWhWXZs1I2FhxF/oLmiKbY75+UcK6OwFPAAdEZK8p+z/gU2C+iIwA/gEeBFBKHRKR+cBhDA+0F5RSmb92zwGzAFdghXmAYcDmiMgJjJHLw2ZdUSLyAbDDvO99pZR18p0WkfV/X+bNBfuIjE/hjV6NeaZLvSLta7E2wYHeLHupMx8vP8K0jafZeOIKEx9uTeManrZWTVMBOBVzive2vEeb6m14qc1L+d6XuGsXKWfOEDByZBlqVzZ4OnnySNNHmLJ/CiejT1Lfp77F6nYPCSFy2jTS4+Kw96x4/5P5/sIppXaZbsJPK6XW5z4Kq1gptVEpJUqplkqpVuaxXCkVqZTqoZRqaL5GZSvzkVKqvlKqsVJqRTb5TqVUsHntRXOUglIqSSn1oFKqgVKqg1LqVLYyM0x5A6XUzBL/hUpJYko67/x2kKEztuPl4sjiFzrxwp0NyoVxycTVyZ4PBgQzY1g7rsQn0++/G5mx8TQZGTqemSZ/EtMSeS30NVzsXfi8y+c42uUfQSp6wULs3N3x6t0r33sqMo83fRxXB1emHbDsKMajcwikp5OwZYtF6y0rCvyVM0cQfiJS/BVgDXvPRdN30ga+33KWESF1WfpSCMGBxd80WVZ0b+LPytFd6NKwGu//fpihM7dzKVbHM9PciFKKD7d+yMnok3za+VP83fPf05IeF0fsqlV49e2LnZtbGWpZdlRxqcKDjR5kxekVnIs7V3iBIuJ6663YubtX2GmyojxGnwE2icjbIvJq5mFlvSo0qekZfLXmbx6YvJmkVGPT5Nv3lm7TZFlRzcOZqUPa8dHAYHaciaLXhL9YeTDc1mppyhmLTyxmycklPHPrM3QM7FjgvbHLlqMSEyvF3peCGNp8KHZix8yDlpswEUdH3O6ouFkui2JgLmBEUbYDPLMdmjw4ERHPA5M3M+nP4/S3wqbJskBEeOy2W1j2cmdqV3Xj2R928+aCfcQn63hmGjgWdYyPtn3EbQG38WzLZwu9P3rhQpwbNcLFyq62tqa6W3UGNBjA4hOLuZRwyWL1eoSEVNgsl4Xug1FKvQcgIu5KqQTrq1SxWLznPONWHeN8dCLe61cRn5SGl6sjkx9rQ58WFTttTn0/DxY+15GJfxznm9ATbD0VxfiHWnEu6lpWnwO3ruWNXo3zjOysqXzEp8Tz2vrX8HLy4tPOn2JvV/CoPOnY3yQdOID//42pkIEti8vw4OH8evxXZh+ezZvt37RIndezXG6ocB54RQl2eYeIHMYI9YKI3Coi31hdswrA4j3nGfPrgawUxjGJaSjglbsaVnjjkomjvR2v92rMvJF3kJ6hGDR5M6//si+rz+ejExnz6wEW7zlvY0011kYpxdgtYzkXd47Pu3xONdfCR+bRCxcgjo549etXBhraniDPIO6pew8L/l7A1aSrFqnTKSjIyHJZAcPGFGWKbAJGbK9IAKXUPqCLFXWqMIxbdYzEXPtGMhR899dpG2lkPTrUrcqK0Z1xcbQnLZd3WWJqOuNWHbORZpqyYt6xeaw6s4qXWr9EuxrtCr0/IyWF2N+W4HFXDxyqlOtYsxblqRZPkZSWxJzDcwq/uYhkZblMTrZYnWVBkXxllVK53SL0bjzggvkUX1R5RcfLxTHfjZiVtc8ag4NXDvL5js/pEtSF4cHDi1Qm/o8/SI+JqXBZK0tLPZ963HXLXcw7Oo+4lDiL1OneOQSVlFThslwWxcCcE5GOgBIRJxF5HXO67Ganpo9rseSVgZuxzzc7MckxvL7+dfxc/fio00dFzuAYvWAhjjVrVvisjCXhqRZPEZcax8/HfrZIfe4dOiCOjsRXMHflonxTngVewIhGHAa0Mt/f9LzRqzGuuVyPXR3teaNXYxtpZH3y6jNA10YVy1NOUzDLTi2j54KevHT2JXr80oML8Rf4ousX+Lj4FKl8Sth5EjZvxvv++xErRwIujzTzbUanwE7MOTyHxLTSj+7t3Nxwbdu2woXvL0qwyytKqceUUv5KqepKqceVUpFloVx5Z0DrQD65vwWB5tN7oI8rn9zfolJ7VOXuc4C3Cw2ru/PT9nN8/efxCumrr8nJslPLGLt5LOEJxv6n5PRk7O3si7WBMObXX0EEn/sHWkvNcs/IFiOJSopi4d8LLVKfRwXMclmom7KIzCTvUPdFm4it5AxoHciA1oE3VSj33H1OScvgrYX7+XLN35yPTuSDAcE4lqNQOJriMXH3RJLSc0ZwSMtIY+LuifSt17fQ8io9nehFi3Dv1AnHmjWtpWa5p41/G9r6t2XmoZkMbjwYJ/vSBURxDwmBL74kYeNGfB6430JaWpei/Ar8Diwzjz8BLyDemkppKhZODnZ8OfhWXryzAfN2nOOp2TtJ0JsyKywXEy4WS56bhM1bSAsPr/Q794vCyBYjibgWwZKThQagL5SsLJebKs40WVGmyBZmO34EBgPB1ldNU5EQEV7v1ZiPB7Zg44krPDRlCxFxOo5ZRaSGe41iyXMTvWAB9lWq4NG9uyXVqpDcUfMOmvs2Z/qB6aRllO6hS0Tw6NiJhE2bUekVw5G3JPMYDYHallZEUzl49LbaTBvSjpMRCQz832ZORFjGTVNTdjTwvjH1r4u9C6PajCq0bFpUFHFr1+J9333YOekYuSLC0y2fJiw+jFVnSp+Syr1zZ9JjYkg6ZNnkZtaiKDv540QkNvMVWAr8y/qqaSoqdzapzs/P3E5yWjoPTN7C9tM2TcWjKQarz6xmw4UN3F7jdgLcjWgUAe4BjO04tkjrLzFLlkBqqp4ey8adte6kgU8Dph2YRobKKFVdWVkuK4g3WVGmyDyVUl7ZXhsppSzjFqGptLQM8mHR853w9XDi8Wnb+H1/uclYrcmHUzGneHvT27Ss1pL/3fU/Vg9azde3fM3qQauLtrivFNELFuBya0ucGzYsA40rBnZix4gWIzgRfYJ159aVqi6HKlVwad68woTvL8oIpk1BR1koqamY1Krqxq/PdeTWWt68+NMepv51Srsxl1OupV7jlXWv4GzvzJfdviyRx1PSvn2knDh50+3cLwq96/QmyCOIqfunlvp/wD2kE4n79pEeG2sh7axHUdZgvgG2AlOAqcA2YBLwJSVPp6y5SfBxc2LOiNvo2yKAj5Yf4b2lh0nXmTLLFUop3t38Lmdiz/B518+LvJifm+iFCxE3N7z63GNhDSs+DnYOjGgxgkORh9hyoXTZKT1CMrNcbrWQdtajqAnH2iql2iml2gKtgRNKqTuVUtpNRFMoLo72fP1Ia54KqcuszWd4/sddJKZUDC+Ym4Efj/zIyjMrean1S9wecHuJ6shISCB22XK8evfG3sPdwhpWDu6rfx/V3aoz5cCUUtVzPctl+V+HKYqBaaKUOpD5Ril1ECNcTIGIyAwRiRCRg9lkY0XkvIjsNY97sl0bIyInROSYiPTKJm8rIgfMa5PETCohIs4i8rMp3yYidbKVGSoix81jaBH6qLEydnbCf+5txrv9mrH68CUenbaVyPiKFRm2MrInYg9f7vySbrW6FTmIZV7ErlxJxrVrenqsAJzsnRjWfBi7Lu1i96XdJa4nM8tl/Kbyn+WyKAbmiIhME5FuItJVRKZStGCXs4DeecjHK6VamcdyABFpBjwMNDfLfCMimQGvJgMjMdyjG2arcwRwVSnVABgPfGbWVRV4F7gN6AC8KyI3T6zwcs6Tneoy+bE2HL4QywOTN3Pmis5hZyuuJF7htdDXCPAI4KOQogexzIvoBQtxqlcP19atLKdgJeSBhg9QxbkKUw9MLVU9HiGdSbsQTsrp8p0apCjfqCeBQ8AoYDRw2JQViFLqL6Co/qn9gXlKqWSl1GngBNBBRAIAL6XUFmWY6u+BAdnKzDbPFwA9zNFNL2CNUipKKXUVWEPehk5jI3oHB/DT07cRk5jK/ZM3s+cfyyRm0hSdtIw03lj/BnEpcYzvNh4vJ68S15V88iSJe/bg88ADN0XWytLg5ujGE82eYOP5jRyOPFzieq5nuSzf02RFcVNOUkqNV0oNNI/xSqnSbNF+UUT2m1NomSOLQCB7JL0wU5YZwTm3PEcZpVQaEAP4FlCXphzR9paqLHyuIx7ODjwydSurDxUtDInGMkzcPZGdl3byzh3v0Lhq6aJ/Ry9YCA4OeA/obyHtKjcPN3kYT0dPph2YVuI6nIICcapTp9zvhyk02KWFmQx8gBE88wMMT7ThQF6PPaoAOSUskwMRGYkx/Ya/vz+hoaEFqF4w8fHxpSpfEbFEn1+/Fcbvhmfm7OKxpk7cdYujZZSzApXlM96bsJdZV2bR2aMznuc8CT0XesM91S+tp96pOXRNvkzSFj9O1XuCCP+uN1aWlobfgl9IaRHMxgMHbrxeASmLz7mja0dWn13NvDXzqOFYMq89z7p1cN24idA1a8CxdP831upzmRoYpdSlzHNzLed3820YUCvbrUHABVMelIc8e5kwEXEAvDGm5MKAbrnKhOajzxQM92vatWunShMN+WaKppyJpfp8151pvDx3Dz8cicC9ehD/6tUEO7vyN9Vik894/3z4832ICQPvIOjxDrQcXOLqTsec5q1lb9GyWksm9J6Q936X/fNh02RINfKYuCRfptmJyTRr2vSGtmNXr+Z8XDwNn30Wj655GKAKSFl8zrcm3cpfC/9iv8t+Hu78cInqiAPC1oXSzs0Nj06dSqWPtfpcZAMjIl6AUkqVOLiUiAQopcLNtwOBTA+zJcBPIvIVUBNjMX+7UirdDFFzO8b+myHA19nKDAW2AIOAtUopJSKrgI+zTb/1BMaUVGeN9XFzcuC7J9rx7pKDfLf+FBeik/jiwZY4O9yY2OymYv98WPpy1g89MeeM92D80KenQloSpCZBWqL5ah6pidnOjevXUuJ45fR8nNLT+NK+Jk6r38lWLhHSko1y/2yB9JScuqQmwsq3ILAtVKkDdsZnE71wIQ7+/llrApqiUcWlCg80fIC5R+fyfKvnCfIMKrxQLjKzXCZs3FRqA2MtipIPph0wE/A03ko0MFwpVWByaBGZizGSqCYiYRieXd1EpBXGlNUZ4BkApdQhEZmP4UCQBryglMrcKPEchkeaK7DCPACmA3NE5ATGyOVhs64oEfkA2GHe975SSgfDKufY2wkf9A8m0MeNz1YeJSI2iSlPtMPbrfxOmVmUjAxIiICY84YhiT0P6z6+blwySU2EX0fComdBFX0vkQLG+vly2t2N7y5GUOP0THB0BQcXcHQxXh1cDFlu45LJtUj4ug04uIJfY1Kd65Hw1zZ8B/dE4i+CV03Qi/xFZljzYfx87GdmHpzJ23e8XezyObJc/utNK2hYeooygpkBPK+U2gAgIiEYBqdlQYWUUo/kIZ5ewP0fAR/lId9JHukBTEeDB/Opa4apt6YCISI8160+NX1ceP2XfQz6djMzn2xPUBU32ypmTlN1jQmDPSWcpkqKNaa5YsIg1nyNOZ/t/XnISC1iZQpCRhs/9LmNg4NznvKfzq5kxYFvGdXyOW4f8kzWCCRPxgcbRi43HtWhx7sQcQQuHSJm1WZQdvgkfA/jZ4KLN1RvBtWbmq/muVvV4v2tbhL83f3p36A/i04s4plbn6G6W/Vi1+HROYSIcV+QeukSjv7+VtCydBTFwMRlGhcApdRGEdEx2DVWo3+rQKp7ujByzk4GfrOZmcPaExzobRtlsk1TCdw4TQWQlmKMOGJNgxFzLpvxMF+Tc8WNEnvwCgTvQAhqD80GGGss2Y9vO+f9Q+9dyzByRWRPxB6+ODjN2EzZ6lkobL9Lj3dyTs2BYbx6fpTVZ5WRQfT3vXBr74/Ti89BxGHzOAIHFkJytuc7jxrXjY6/aXT8moBTPjv+LbzuVJ4ZHjycRccXMfvQbN5o/0axy7uHhMC4L0jYuKlcZrnM18BkC2S5XUS+A+ZijLQfIp9Fc43GUtxR35eFz3Vk2IztPPTdFr55vC1dG/mVvSJ/vp/3NNWSl2Hbt8aPYHwENzgqulUzjEfVelCns2k0Ag3j4B0EHv4FjyIg/x/6YhiXEm2mzPwx//N9VEwYkseP/LXt20kNC8Nv1Cio08k4MlEK4sLhUjajE3EYdk431oQAEKhyC1Rvbhof0wCF74Vlr+a/7lTJqOVZiz51+zD3yFxWnVlFxLUIarjXYFSbUUWKYO3cqBEOfn4kbCqfaZQLGsF8mev9u9nOy3d8Ak2loJG/J4te6MSwmTsYPmsHD7YNYsPxK1yITqSmjytv9GrMgNZW2uKUkQEX9+U9ggBjUTxzSsi7lmk8goxzr5qGISgt2X7oS/I0n30z5eS7JhdvM2XLwdByMOvz8S6K/mUBdt7eePa8+8ayIsbfwKsmNLzrujwjHa6eyWl0Lh2Gv1cWvJ6Ummj8DSqhgQFo6NOQ39XvXLpmONmGJ4QzdvNYgEKNjIjg3qkT8evWodLTEfvy5RiTr4FRSt1ZlopoNHnh7+XC/GduZ9Dkzczbcf3H/nx0ImN+NfZdWMzIJMfByXVwfBUcXwPxl/K/17sWPLHIMu0WhPlDXxIyN1N+HPJxqTdTZic9Joa4NWvwefBB7Jydi17Qzh586xtH037X5WnJcOW4YXR+fSrvsjFhecsrAfOOzbtBlpSexMTdE4s0inEPCSFm8WKSDh7E9dZbraFiiSmKF5kz8ABQJ/v9Sqn3raeWRnMdTxdH4pJvzGeemJrOuFXHSm5glILIE/D3Kji+Gs5uNhbanb2hQXdo2AtSrsGaf5dqmsoWrDm7hlmHZvFQ44foV79f4QWKQczS31EpKZbLWungDDWCjePP9/IeNbr6GKNKu5LHSyuvXEzIO4pFfvLcZM9yWeEMDPAbRhiWXYAOf6uxCeHReUcnuhCdmKc8X9KS4cxGw6D8vQqumsEC/ZrCHc9Dw55Q6zawz+Ye7eJZ4HpEeeN0zOmszJRvtres+2pW1srmzXFp2tSidQN5rzuJHSRehVl94d7xUL2J5du1ITXcaxCeEJ6nvChkz3Lp98ILllavVBTFwAQppXSwSI1Nqenjyvk8jImbkz0JyWm4OxfwVY69YBqU1XAqFFITDPfdul3gjhcMo1LllvzLF7IeUZ7IzEzpZOdU4syUBZF06DDJR49S410rjeDyWnfq/jakJ8Pqt+HbEOg0Crq8bpl1rnLAqDajGLt5LEnp1x+iXOxdGNVmVJHrcO8cQuSUqaTHxmLvVfLApZamKAZms4i0yJ4TRqMpa97o1Zgxvx4gMfX6YrC9nZCQkk6fiRv4cvCttK9j7rfISIewncZayt+r4ZL51fWuBbc+DI16GZ5dTjbeX2NhlFKM3TyW07Gn+e7u70qcmbIgohcuQJyd8epb+NpAiclv3alRH1j9H9jwBRxcCPd+BfUrfs7DzHWWibsnEp4Qjr3Y887t7xRp/SUTj5AQIid/S8KWrXj16mktVYtNUQxMCDBMRE5jTJEJRsiYAjdaajSWZEDrQALP/U6t3eOori4TIX6ca/MGGcGDeH3BPp7+bg1jm12kn+sB7E/+CYlRxl6T2rfDXe8Zo5TqTSv1TvOfjv7EijMrGNVmVIkzUxZERmIisb8vw7NXT9s8JXv4wf3fQatH4PdXYc5AaPEg9PrY2ARagelbry996/Xlj7N/8EroKzg6FC+ChWvLlth5eJCwcWOFMzB9rK6FRlMY++fT/sC7QCII1OAyNfa/DUlbCK16Abm2A7uTGUSLF6r+XVRpda/xdOt6c+Sa2xOxhy92fFHqzJQFEbd6NRlxcbbPWlmvGzy3GTZ+BRvHG9Ofd70HbYZWeCeA7rW7U8erDjMOzKDXLb2KnF9HHB1xz5blsrzk5SlKPpizeR1loZxGk0VeGx7TkuHwIuzTk7Dr8jq7755PH4dptD/8IBMutiDVyUa7/8uYK4lXeD30dYtkpiyI6AULcbylNm7t21ul/mLh6AJ3/h88uwn8W8Dvo2Fmb7h0yNaalQo7sePJ4Cc5EnWELRe2FKuse6cQI8vlqVNW0q74VGxzr7l5yHcfhMAzf0H3f9OmUy9WvtKde1sGMOGP49z/zWaOX6rcUY3SMtJ48683iU2JLXVmyoJIOXOGazt24HN/Octa6dcIhv0OAyYbe2m+6wJ/jDXcyyso99a7l+qu1ZlxsHjhFMtjlkttYDTlm8vH4OcnyDd4hHfOMOfebo5MeLg1kx9rw/noRPp+vZGpf50iPaNyBp+YtHsSOy7usEhmyoKIXvgr2NvjPWCA1dooMSLQ6lF4cSe0fNiYNvvmdmOzbAXEyd6JIc2HsO3iNg5cLrpv1fUsl5usqF3x0AZGUz6JPgeLXzB+KE6uhSb9jCjB2Slgw2OfFgGsGt2Fro38+Gj5ER6ZspV/IivuU21e/HH2D2YemmmVzZTZUWlpxCxejEeXLjj6l+PFdHdfGPA/GLbM2Lz54yD4ZRjEVbx03IMaDcLTybP4o5jOnbm2YwcZyeVjy6I2MJryRcIVWDnGyDty4Be47TkYtQ8e/gHum2S4GiPGa79JBW549PN0ZsoTbfnywVs5Eh5L74l/8eO2syhV8Uczp2NO859N/7HKZsrcxP+1gbTLly23c9/a1AmBZzfCnf+Bo8vhv+1h+1TDfb2C4O7oziNNHuHPf/7kVEzR11Q8QjqhkpK4tnOnFbUrOtrAaMoHSbFGgq2JtxpRils+BC/vht4fg3s1456Wg+GVgzA22ngtwm56EeGBtkGseqULbWpX4d+LDjJ05g4uxuQdGaAiYO3NlLmJXrgQ+2rV8OjSxartWBQHZ+j6Bjy/BQLbwPLXYfrdEL7f1poVmceaPoazvTOzDs4qchm39u2zslyWB7SB0diW1ETY/LVhWNZ/Bg16wPPboP9/b1hfKQ01fVz5fngHPujfnB2no+g5fj2L95yvcKOZ7JspP+/6uVU2U2bHLiaG+NBQfAYOQBwrYHZR3/rwxGK4fypE/wNTuhmbNVMSbK1ZoVR1qcqABgNYemoplxIKCLyaDTs3N1zbtS03C/3awGhsQ3oa7JoFk9oY//A1W8HIUBj8veEZZAXs7IQn7qjDilGdaejvyeif9/L8j7uJjC8f89VFIXMz5UutX7LKZspMYpYu5Xj3HlT711uQno6dTwXeTyRijHZf2A6tHzceaP53GxxbaWvNCmVo86EopZhzeE6Ry3iEhJB8/DipF22/9qQNjKZsyciAg7/C/zrA0lFGHpWhvxuh72u2LhMV6lRzZ/4zd/BWnyb8eSSCXhP+YvUh2/8zFkZZbKYEw7iEv/0OaRcukOmQfOXrr4lZutRqbZYJblWNdbzhq8DJA+Y+BD8/bsSqK6cEeQbRu25vfvn7F2KSY4pUJstdeZPtp8m0gdGUDUrB8T9gSldY8CTYO8HDP8GINVC3c5mrY28nPNu1PktfCsHfy4WRc3bx6vy9xCSmlrkuBbHs1DJ6LuhJy9ktGbZyGF5OXlbdTAkQMX4CKinnGpVKSiJi/ASrtVmm1L7d2DvV4x3Dlfm/HWDrt7BvHowPpmvoABgfbKRuLgcMDx7OtbRrzDt6Y96YvMjMchlfDqbJrPYtFZEZIhIhIgezyaqKyBoROW6+Vsl2bYyInBCRYyLSK5u8rYgcMK9NEnOXl4g4i8jPpnybiNTJVmao2cZxERlqrT5qisg/24xQ6z8+AEnRMPA7eG4TNOlr89hgjWt4suj5TrzcvQG/7b1A7wl/sfH4FZvqlMmyU8sYu3ks4QnhKBQZKoP4tHg2hG2wartp4TeGji9IXiFxcILOr8HzW6FWB1j5L1j8LMScQ1DXUzWXAyPTqEojugR14ccjP5KYVnh6ChHBPSSEhM1bUOm29Zyz5ghmFpA7zP9bwJ9KqYbAn+Z7RKQZ8DDQ3CzzjYhk5v6cDIwEGppHZp0jgKtKqQbAeOAzs66qGOmdbwM6AO9mN2SaMuTiQfjpIZjR09hlfc8X8OIuI6JxYfnoyxAnBzte7dmYX5/riJuTPY9P38Y7vx3kWsqNSc7Kkom7J+YI4Q6Qkp7CxN0TrdquQ0BAseQVmqp14fGF4OZrjLKzk5mquRwwIngEV5Ovsuh40bKouod0IiMmhqSDBwu/2YpYzcAopf4ConKJ+wOzzfPZwIBs8nlKqWSl1GngBNBBRAIAL6XUFmW4+3yfq0xmXQuAHubophewRikVpZS6CqzhRkOnsSZRp2Hh00bujrNbjKmIUXuhw9PGk2M55dZaPix7uTNPhdRlztaz3DNxA7vO5v4Klx2lzXRYUjx69LhBJi4uVH9ltFXbtRkicC2fz7mcpGpu49+GVn6tmH1oNqkZhU/june8nuXSlhQlmrIl8VdKhQMopcJFJHNbcCCwNdt9YaYs1TzPLc8sc86sK01EYgDf7PI8yuRAREZijI7w9/cnNDS0xB2Lj48vVfmKSO4+OyVHccvZ+QSEr0aJPedrDeSf2veTlu4Jm3fYTtFiEuIB1dq7MO1AIoMmb6FlNTvOxSuikjLwDV3OA40c6VjTui67qSoVB3EgVd34Y+Jj72O175rdlUh8Fywgo1o1JD0du6tXyahalfj+/bno6QmV9Dt+u3M1XJIv3yBPcq7G1nLS5w50YErCFCYsn0B7j8IDjlatXZsLy5ZzqHnzQu+11u9XWRuY/MhrIl4VIC9pmZxCpaYAUwDatWunSpOtMLQCZDu0GPvn50wh3Pk1uHoGdnxn5LRvOwy6vkltzxrUtrWuJaQb8HjfNJ6evYMtpzKfboXIJMWcI+k0a9qMAa3zfG4pNbEpsYxaO8owMnYOpGVcn6pzsXfhXx3/Rbd63SzerkpN5ezjT5Bsb0/deXNxCgq6eb7XVT++MVUz4FI/pNz0v4vqwtola9mcvpnXu75eaNDRiH37iPxuCp1bt8beu+DI4tb6nMvai+ySOe2F+RphysOAWtnuCwIumPKgPOQ5yoiIA+CNMSWXX10aS7B/vvGPmH0x9PfRsGkCNO0HL+4wMg16WncDYFng4ezAP1E3LqompqYzbtUxq7R5KeESw1YOY+/lvXzW+TM+7PQhAe4BCEKAewBjO44tVqbD4nB50iQS9+0j4MMPcAqy3CbXCkHLwUboIe9aKMTY5BvUAQ4vgm3f2Vo7wAjlPzx4OCeiT7DhfOGOHh4hIZCRQcKWrYXeay3K2sAsATK9uoYCv2WTP2x6htXFWMzfbk6nxYnI7eb6ypBcZTLrGgSsNddpVgE9RaSKubjf05RpLEFeeVkAPPzhgalQtV7Z62RFLkTn7bVzPjqRnWcsuz5zKvoUT6x4gvNx5/mmxzfcU+8e+tbry+pBq9k/dD+rB622mnGJ37iJyKnT8Bk8GK/eN+mSpRmKaH23xfDKIXhyOTS5F1a8CbuLvtHRmvSu25sA9wCmH5he6L1ZWS432W4dxppuynOBLUBjEQkTkRHAp8DdInIcuNt8j1LqEDAfOAysBF5QSmX61z0HTMNY+D8JrDDl0wFfETkBvIrpkaaUigI+AHaYx/umTGMJ8lv0jI/IW17BqenjmqfcTmDQt1t48NvNrDsaUeqQM3sj9jJk5RBS0lOY1XsWd9S8o1T1FYe0y5e58K9/4dywIf7/N6bM2i332DvCoBlGZtQlL8GBBbbWCEc7R4Y2H8ruiN3sidhT4L1Glss7iN+4yWYhkazpRfaIUipAKeWolApSSk1XSkUqpXoopRqar1HZ7v9IKVVfKdVYKbUim3ynUirYvPaiOUpBKZWklHpQKdVAKdVBKXUqW5kZpryBUmqmtfp403FsRf77ViwYN6w88Uavxrg65nSpdnW059P7W/Buv2acv5rIk7N20GfiBn7be5609IxitxF6LpSnVz+Nj7MPc+6ZQ1PfphbSvnBURgYX/vUvMhISCPzqS+xcXMqs7QqBgzM89CPUvgMWPWNEZ7YxAxsMxMfZhxkHCg/l7x4SQlq47bJc6p38msJJS4YV/4K5D4NnIDjk+hEqIC9LRWdA60A+ub8FgeZIJtDHlU/ub8Hg9rV5slNd1r95J18+eCvpGYpR8/Zy55ehzNl6lqTUom1wW/j3QkatG0UDnwZ83+d7annWKryQBYmcOo2EzVvw//f/4dywYZm2XWFwcoNHf4YaLeGXoXBynU3VcXN049GmjxIaFsrxq8cLvNcjpBNguyyX2sBoCubKcZjWwwihf9tz8PIuuO/rbIuhhedlqegMaB3Ipre6M6u3O5ve6p7De8zR3s5IBzC6C1OeaIuvuzNvLz5IyGdr+d+6E/mGnlFK8e2+bxm7ZSwda3Zkeq/pVHWpWlZdAuDa7j1cnjQJr3vuwWfQoDJtu8Lh4mVsyPRtCPMeNfZ32ZBHGj+Cq4MrMw8WPEHjGBiIU926NstyqQ2MJm+Ugr0/wXddIeY8PPIz9PnUmDLIsRhatLwslR07O6Fn8xoser4j80beTvOa3oxbdYyQT9fy6YqjRMRd35GfnpHOh1s/5H97/8d99e9jUvdJuDm6lam+6dHRnH/9NRxr1qTGe2MLdXnVYATLHLIYvGrCT4PhQsFrINbEx8WHQY0Gsfz0ci7EF+wk6x4SwrXt28nIFV+uLNAGRnMjSbHw60hY/JyRrOm5TdD4JvUsKiYiwu31fJk9vAO/vxRC18Z+TPnrJCGfreP/Fh3g74hIXg19lfl/z+epFk/xYacPcbQr2zwrSinC336btMtXCPzqS+w9Pcu0/QqNR3UY8hu4+MCcgXDpsM1UGdJsCCLC94e/L/A+j5BOqORkru3cVUaaXUcbGE1Ozu+G77rAwQVGytkhvxlPbJpiExzozX8fbcPa17oxqG0QC3b/zYCFQ1l7bh3DGo9mVJtRNhk5XJ07l7g1f1D91VdxbdGizNuv8HgHwdDfwN4Z5gyAyJM2UaOGew3urXcvC/9eyNWkq/nedz3LZdmvw2gDozHIyDASMU3vCempMGy5kXK2HAWlrKjUqebOyz39aNR6No5u51GXHufrxTUYOmM7W09FlqkLadKRI0R8+hnuXbtQdeiQMmu30lG1nvHwlZEGs+8zsmXagCebP0lSehI/Hf0p33vs3Nxwa9/OJvthtIHRQPxl+OlBI7Nko17w7Aa4pez2YVR2jl89zmPLHyMy6TLTek1h00uv8Eavxhy6EMPDU7bywOTNrDl8iYwM6xqajIQEzr/6GvY+PtT85BPETv/7l4rqTYxEeclx8H1/iCv7pHX1fOrRvVZ3fjryE9dSr+V7n3unEJKPnyjzLJf6G3azcyoUvu0EpzdA3y/hoR+MxUyNRdh5cSdDVxppb2f1nkX7Gu3xdnXkhTsbsPFf3fmgf3Mi4pJ5+vud9JrwFwt3hZFagr00ReHihx+RcuYMNceNw6Gq/owtQsCt8PgCiLsE3w+AhMgyV2F4i+HEpsSy4O/8N4LaKsulNjA3K+mp8Md7xj+Fiw+MXAftn7J5ArDKxB9n/+CZNc/g6+LLD/f8QOOqjXNcd3G054k76hD6ejcmPtwKezvhtV/20W1cKLM2nSYxJZ3Fe87T6dO11H1rGZ0+XcviPedLpEvMkiXELFpEteeew/22DpboniaTWh3g0Xlw9TT8MBCSipba2FLc6ncr7Wu0Z/bh2aSm5+0W79yoIQ7Vq5d5+H5tYG5Grp6FmX1g41fQZgiMDAX/wkN6a4rOz0d/5tXQV2nq25Q5feZQ0yN/RwkHezv6twpkxajOzBjWjpo+Loxdeph2H67h9V/2cT46EYUR/2zMrweKbWSST58mfOx7uLZrS7XnnytlzzR5UrcLDJ5jeJX9OBhSEsq0+RHBI4i4FsGy08vyvC4iuHfqVOZZLrWBudk4tAi+7QyXj8GgmXDfJGOnssYiKKX4es/XfLjtQ7oGdWVqz6n4uPgUqayI0L2JP78825Ffnr2DtAxFWq51meJGcs5ISeH8a69h5+hI4BdfIA7lJUNHJaRRT3hgGoRtNzZjppbdvpOONTvSpGoTZhycQYbKe4o1K8vlgQNlppc2MDcLKddgycvwyzCo1tBYyA++39ZaVSrSMtJ4d/O7TNk/hQcaPsD4O8fj6pB3sMzCaF+nKilpef9QnI9O5M8jl0hOK/xJNGLcFyQfPkLAJ5/gWKPip1Ao9zQfAP3/Z6xt/jLMmIouA0SE4cHDOR1zmnXn8g5lcz3LZdmtw2gDczNw6TBMvRN2fw8hr8DwlVCljq21qlQkpiUyet1oFp1YxLO3Psu7d7yLg13pRgv5RXIWYMTsnbT94A9GzdvDyoPhJKbcaGzi1q7l6pw5VBnyBJ7d7yyVLppi0OpRuOcL+HuFsWE5o2ympO6+5W6CPIKYfmB6nq7vDlWq4BIcXKb7YfR4uTKjFOycAav+D5y94IlfjdDjGotyNekqL659kYNXDvL27W8zuLFlQue80asxY349QGK2wJmujvZ80L851TydWXHgIqsPX+S3vRdwdbTnziZ+9AkO4M4m1XGOukz4mP/DpVkzqr/+ukX00RSDDk9D6jVY8w44uhnx+6zsFu5g58CTwU/ywdYP2HlpJ+1r3JhW2aNzCFe+/Y70mJhCs1xaRCert6CxDYlXjRwWR5ZC/R4w8FsjzIXGopyPP8+za57lQvwFvur6FT1u6WGxujODao5bdYwL0YnU9HHljV6Ns+TdGlfno/Rgtp2OYsXBcFYevMTyAxdxsVN8vX0qAckp+H78GXZOThbTSVMMOo0yFvvXfwZO7tDnM6t7afZv0J9v9n7D9IPT8zQw7iEhXPlmMglbtuLVu5dVdQFtYCon/2yFhU9BXDj0/BBuf8HqT083I8eijvHcH8+RlJ7E1J5TaePfxuJtDGgdmCN6c24c7O3o1KAanRpU4737gtl19ir/fDGeoLC/+bzto2yc8zedGkRxT3AAdzfzp4q7NjZlSrcxhpHZ8l/DyNz1rlWbc7Z35vFmjzNx90SORh2lSdUmOa67tmyJnacnCZs2lomB0b86lYmMdFg/DmbeY4R4GbEaOr6kjYsV2B6+nWErh2Endnzf+3urGJfiYm8nNL90nOZrF+I1cCAvffIiT3aqy8nL8by5cD/tPvqDx6dt44etZ7kcl2xrdW8ORIyHvLZPGtsC/vrC6k0+1PghPBw98kxIJg4OuN9+O/EbNpZJiCI9gqksxIbDr0/DmQ3Q4kHo+5WRw0JjcVaeWcn/bfg/anvW5tu7v6WGe/nwzkqLjOTCG2/gVKcOAW//h0A3N9rUrsKYPk04dCGWFQfDWXHgIv9ZfJC3fztI+zpVuSe4Br2DA6jhrTNZWg0R4/8xJQHWfgBOHnD7s1ZrztPJkwcbP8jsQ7N5KfYlannlTGLnHhJC3Jo1pJw8iXODBlbTA7SBqbjsnw9/vg8xYeDmC6mJgIL+3xheLHpHvsVYdmoZE3dPJDwhHK+5XsSmxNKmehsmdZ+Et7P1F0qLgsrI4MKYMaTHxFBr2lTs3K7vbRIRggO9CQ705vWejfn7UjzLD4Sz8uBFxi49zNilh2lT24c+wQH0Dq5BrapuLN5znnGrjnE+OpHArWtzrP1oSoCdHQyYbCz8r/yXsfesjfWCjT7R9Al+OPwDsw7N4u073s5xLTPLZfzGjdrAaPJg/3xY+rJpVIBrVwCBu9+D1o/ZVLXKxrJTyxi7eSxJ6camudiUWOzEjgENBpQb4wIQNWs2CX9toMa77+DSuHG+94kIjWt40riGJ6/c3YiTl+NZefAiyw+E89HyI3y0/AhBPq5cjE3K2uSZGUEA0EamNNg7wKAZMPcRY0+aoxu0sE4mUT83P/o36M/iE4t5rtVzVHOtlnUtM8tlwsZN+A4bZpX2M7HJ5LyInBGRAyKyV0R2mrKqIrJGRI6br1Wy3T9GRE6IyDER6ZVN3tas54SITBIzuYaIOIvIz6Z8m4jUKfNOWpM/379uXLJQsH2qTdSpzEzcPTHLuGSSoTKYvG+yjTS6kcT9+4n46is8774bn4cfLlbZ+n4evHBnA5a93Jm/3riTMX2acCkuKc8IAp+vOmpJtW9OHJyNgLK3dDT2yBzNO7SLJXiy+ZOkqTR+PPLjDdfcO4dwbccOq2e5tOXq751KqVZKqXbm+7eAP5VSDYE/zfeISDPgYaA50Bv4RkQyk5RMBkYCDc0jM+3iCOCqUqoBMB74rAz6UzZkpEPMubyvxYSVrS6VnDMxZwhPCM/z2sWEsg/NnhfpcXGcf/U1HKtXJ+DDD0qVwKy2rxvPdK1PWnrei78XopN4df5elh8IJy6pbHaoV0qc3ODRn6FmK2O3/8m1Vmmmtldt7r7lbuYdnUdcSlyOax4hIWWS5bI8uRf1B2ab57OBAdnk85RSyUqp08AJoIOIBABeSqktynCH+D5Xmcy6FgA9pDT/eeWFyJNGkMr88A4qO10qMVFJUXy87WMG/jYQIe+vTXlY2FdKEf7OO6SGh1Pzyy8stnEuvwgCro72rD0awfM/7qb1+2t4bNpWpm88zZkrZRvYsVLg7AmPLYBqjWDuo7D2QxgfDGN9jNf98y3SzPDg4cSnxvPL37/kkLu1b484OVl9V7+t1mAUsFpEFPCdUmoK4K+UCgdQSoWLSOauwEBga7ayYaYs1TzPLc8sc86sK01EYgBf4Ep2JURkJMYICH9/f0JDQ0vcofj4+FKVLxClqHlhBfVPziLDzoFLNfsQcHEt9hnXXU3T7Zw5VvNBIqylQx5Ytc82ICUjhdC4UNbErCFFpdDRoyMBjgEsjl5Mqrr+xO4ojtztcrfN++66YSNeK1YSN2AA22NiwEL69K2dzqxYSMkWCs3JDp5oas9tNRw4GWPP3oh09l2M4oMTkXzw+2FquAut/BxoVd2eBj52ONhV3Oe5svxeO9Z/g3ZRo3D6a9z1R5mYc6QvfpFjR44Q4d+11G00cWnC9L3TqXW5Fo7imCX3qVePiFWrOHL7bVbrs60MTCel1AXTiKwRkYImd/P6pqoC5AWVySkwDNsUgHbt2qlu3boVqHRBhIaGUpry+RJzHpa8aAyj63fH/r7/EuQdmNOLzDsI+x7v0KzlYJpZXoN8sVqfy5gMlcGyU8uYtGcSFxMu0q1WN15p+wr1vOsB0PZU2ywvsgD3AEa1GUXfen1tqnPy8eOcXrAAt44dafLxRxbNTtkNaJbdiyxXBIHssQrORV3jzyOX+PNoBGtPRbHyTCpeLg50bVydHk2q07WRX4Xb3Fnm3+sDb0FqdA6RfUYyzS78QrOHSr8x0zXcladWP0V0YDQPNnowSx558hQR48bRqXFjNh07ZpU+28TAKKUumK8RIrII6ABcEpEAc/QSAESYt4cB2R25g4ALpjwoD3n2MmEi4gB4A1HW6o9VUAoO/ALLXzcisvb9EtqNuO5+3HKwcWhKxY6LOxi3YxxHoo7QzLcZH4d8fEOIjb71+tK3Xt9yY1AzEhM5/+qr2Hl6UvPzz6yS+jgzgkBhfa5V1Y1hneoyrFNd4pPT2Hj8CmuPXmLt0css3XcBO4G2t1ShexN/7mpanQbVPUq1TlQpyS/VsoXWVDvU6ECwbzCzDs7i/gb3Y29nLGG7h4TAuHFGlstq1QqppWSUuYEREXfATikVZ573BN4HlgBDgU/N19/MIkuAn0TkK6AmxmL+dqVUuojEicjtwDZgCPB1tjJDgS3AIGCtKottq5Yi4Qr8/gocWQK1bjP8533r21qrSsWp6FOM3zWe0LBQarjX4JPOn3BP3Xuwk/K0LJk3lz7+hOTjJ6g1fRoOVvphKAkezg70Dq5B7+AaZGQo9p+PYa05uvls5VE+W3mUWlVd6dHEn+5NqnNbvao4O9gXXnFlxzsob8cdr/yT1BUHEWFEixG8EvoKf/zzB73qGI6417NcboIB/S3SVm5sMYLxBxaZTzEOwE9KqZUisgOYLyIjgH+ABwGUUodEZD5wGEgDXlBKZYaXfQ6YBbgCK8wDYDowR0ROYIxciue7aUuOLjf2uCTFwF3vmaFe9D+hpYhMjGTyvsks+HsBrg6ujG4zmseaPoaLQ8XYyR67fDnRv/yC79NP49Gpk63VyRc7O6FVLR9a1fLh1Z6NCY9JZN3Ry/x55BJzt//DrM1ncHOyp3PDavRo6s+djavj5+kMkLXJM68An5WSHu/k3NeWSVoyhO+DgFtL3UT32t2p41WH6Qem0/OWnoiIkeUyJIS4P/+E+/qVuo28KHMDo5Q6BdzwF1NKRZJzejf7tY+Aj/KQ7wSC85AnYRqoCkNSLKwcA3t/AP8WMOQ3ncbYgiSmJfLD4R+YfnA6SWlJDG48mGdvfZaqLlVtrVqRSTl3jvB33sW1VSv8Xn7J1uoUiwBvVx69rTaP3labxJR0tpy6wp9HIlh7NIJVhy4BcGstH2p6u7D2aATJZrK1m2KTZ+ZUd7Y1VVo9CrvnwLS7oc+nRiyzUkwt2okdw4OH887md9gSvoWONTsCIK4uZMTEUP35FzhesybVXxmNdz/LGRu9k788cPovWPwCxIZB59eg61vgULEWRssrGSqD30/9zqTdk7h07RLda3VndNvR1PWua2vVioVKSeH8q6+BnR01v/gCcXQsvFA5xdXJnu5N/OnexB+lFEfC47IcBVYcvHE9IjNNdKU1MJD3mmqHkcZmzN9fgbNb4N7x4OxR4ib61uvLf/f8lxkHZtCxZkdili4lZuFCwPCKSrtwgfC33wGwmJEp/xPOlZnURFjxFszuB/aOMHy1MVzWxsUibAvfxkO/P8S/N/4bP1c/ZvaaycTuEyuUcYlZupTj3XtwtOWtJB04gFf//jgFVZ4fWhGhWU0vXurRkMUvdMpn15Exkvlo2WHW/305z+ydlRL3asZemTv/AwcXwNTuEHGkxNU52TsxpPkQtl3cxsErB4kYPwGVlDOqtkpKImL8hFIqfh09grEVYbtg0TMQedx4UrlrrJEvQlNqTkaf5KtdX/FX2F/UdK/JZ50/o3fd3hViAT87MUuXEv72O6hs4TxifvkFt5YtLDqNUZ6o6ePK+ejcYZDAycGO2ZvPMnXDaZzs7Wh7SxVCGlYjpEE1ggO9sa/A+24KxM4Our4BtToYOZ6mdjdGMreWbFl5UKNBfLf/O2YcnMHT4XlHqUjLR14StIEpa9JTYf3nsOFL8KwBTyyG+jpfuiW4kniFb/Z+w8LjC3F3cOfVtq/yaNNHcbZ3trVqJSJi3Bc5jAtcf8KsrAYmvzTRn9zfgp7N/dlx5iobj19m44lIxq06xrhVx/B2daRjfV86NahG54bVqF3VrfK5QtfrCs9ugAUjjAfTs5ugz+fgmHfUhfxwd3TnkSaPMHX/VJ7294OLETfc4xAQYCmttYEpUyKOGF+O8H1w6yPQ+1Nw9bG1VhWexLREvj/0PTMOziAlPYVHmjzCMy2foYpLlcILl0MSDxwgcuo00iJu/OcHyz5hljcKSxPdtZEfXRv5AXA5LpnNJ6+w8fgVNp64krV+E1TFlZAG1QhpWI2O9atRtYJt9MwXzxqG80/ox8YD6vk9MHh2sbcwPNb0Mb4/9D2b+9Wl05zYHA8x4uJC9VdGW0xlbWDKgox02PI/I96Qs6cRTbVp5XwCtSaZeVkuJlykhnsNXmr1Eumk8/Xur4lIjOCu2ncxuu1obvG6xdaqFhulFAmbNhM5dSrXtm3DzssLOw8PMuLjb7jXkk+Y5ZHC0kRn4ufpTP9WgfRvFYhSitNXEth4wjA4y/aHM2/HOUSgeU0vY3TTwI92darg4liB3f7tHYx12lq3w6KR8F1X6P81NB9Y5CqqulRlYMOB/Ff9wt3/eY3kb2aSeuECjtqLrAISdRoWPw//bIYm98K9E8DDz9ZaVThy52UJTwjn35v+jULRoloLxnUdVy7SFhcXlZZG7KpVRE6fTvLhIzhUr071N9/EZ/Bg4tetvWENxtJPmJUFEaGenwf1/DwYckcd0tIz2H8+Jmt0M2Pjab5bfwpnBzva16maNZ3WLMALO3P9pkIlWWvUE57ZAAueNCIyn91ipGYuooPQ0OZDmX9sPr/UucRra/+0WpQKbWCshVKwaxas+rexUXLAt8bCXGWbGy4j8srLolD4OPvw4z0/Vrg594ykJGIWLyZy+gxSz53DqW5dAj76EK9+/bBzMn4kMp8kI8ZPIC08HIeAAIs/YVZWHOztaFO7Cm1qV+HlHg1JSE5j++morBGOEVkAqrg50rF+Ndyd7flt74WKtf/GpxYMWw5/vAtbv4HzO+HBWeBTu9CigR6B9K7bm/nH5vNUi6espqI2MNYgNhyWvAQn1kDdrtD/f8aXQVNsrqVeY0v4lnzzssQkx1Qo45IeG8vVn+YSNWcO6ZGRuLRsSfU338CzR488Y4p59+unDYoFcHd24M4m1bmziRGkPSI2iU0nr7DxeCQbT1zmUmzyDWUqxP4bByfo/QnUvgN+ewG+7QwDv4PGvQstOjx4OMtOLaPPwj7EpcYRsMDygVy1gbE0BxbAsteMMA99Pof2Txuuhpoicz7+POvPreevsL/YfnE7qRmpCIK6MSB2ucjLUhRSL0UQNXs20T//TEZCAu6dO+P71FO4dWhfoQxkZaG6lwsDWwcxsHUQSinqjVmex7fLGMl8tfoYIQ39aF3bB0f7cvq/3Ow+qBEM84fC3Ieg02jo/raxZpMPx68ex07siEs1kpGFJ4QzdvNYAIsZGW1gSosZNr9rzDnY5GpsngxsZzxFVGtga+0qBOkZ6ey/sp/159azPmw9J6JPAFDHqw6PNHmEbrW6cSH+Ah9u/TDHNJmLvQuj2oyyldpFIvnUaSJnTCf2tyWo9HS8+vTB96kRuDRtamvVNCYiku/+G0d74b/rTjBp7Qncney5vZ4vnRtWI6ShH/X93MvXw0HVejBiDax8CzZNgHPbYdD0fINmTtw9kQyVkUOWlJ7ExN0TtYEpF+yfnxWkTsAwLnaOxqhFG5cCiU2JZfP5zawPW8/G8xuJTo7GQRxo69+Wge0G0iWoC3W86+Qo42DnkMOLrDzkZcmPxP37iZw6jbg//kCcnPB5cBBVn3wSp1p6qrQ8UtD+mzubVGfLyStsMB0G/jxquI/X9HYxNns29KNTfV98PcrBfitHF+g3AW7pCEtHG1NmD0zLc69dfmm/LZkOXBuY0vDn+zdGQM1IhXUfQquKE8C5rDgTc4b1YcYoZfel3aSrdHycfegc2JkutbrQqWYnPJ088y2fmZelvKKUImHjJiKnTctyNfZ99hmqPv44Dr6+tlZPUwDZ99/klWStd3AAvYMN9/B/Iq+x4cRlNh6/wsqDF5m/08jbEhzoRUgDPzo3rEbbW2zsDt1ysBGFef4QmDMQur0FXd7IEZm9hnuNPNc2LTntrA1MacgvIZCFEgVVdFIzUtl9aTfrw4z1lLOxZwFoWKUhTwY/SdegrrSo1iIrAVJFJcvVeNp0ko8cwcHfP8vV2N5Dh/+pKBQ1yVptXzce872Fx267hfQMxf6waDYev8KGE1eYtuEU364/iYujHR3q+tLZ3PDZpIZn2U+n+TWGp9fC769C6Cfwz1a4f2rWNolRbUblcP0Hy087awNTGvJLFOQddKOsEpG54TE8IfwGz5OrSVfZeH4j68PWs+n8JuJT43G0c6RDQAcea/oYXYK6EOhRjr1yikFGUhIxixYROWNmNlfjj/Dudy/iVEl2j2sKxN5OaF27Cq1rV+GlHg2JT05j26nIrOm0j5YbwSn9PJ0JMffehDSoRnWvnPmHrJYDx8kdBn4LdTrB8jfgu84waAbc0jHrf9aa6cClIiV6tCbt2rVTO3fuLF6hbGswWTi6Qr9JlTadce4NjwDO9s50C+rGpWuX2Hd5HwpFNddqdA3qSpegLtwecDtujm421Lp0xCxdSsT4CVm7nas9+wzpUVevuxrf2pJqTz+NR/fuVklfbEvKS5rossSSfQ6PSTSMjWlwohJSAGjs70lIQ8PgRMQm8+6SQ3mu/1jURfriAWPK7OpZIxpAx5ezPFxL02cR2aWUapfnNW1gDEpkYCDLi0zFhCHeQcYHV0mNC0DPBT3z3ZPSzLcZXYO60jWoK019m1a46MV5kVdE40zcu5iuxu0rr6uxNjCWIyNDcTg8lo0nrrDh+GV2nLlKSlpGvvcH+riy6a3ullUiKRaWvAiHf4NGfYyIABu+KtXvlzYwRaDEBsaksvwjpmekcznxMmFxYVxIuMD5uPOExYdxPv48F+Iv5GtcBGH/0P1lrK3lyUhIICXsPKlh50j55xxXvv6ajGvXbrjP3s+PRhv+soGGZUtl+V4Xh7Lqc2JKOjvORDFkxvZ87xnQqiYNqntkHbWruuPkUMoHN6Vg+xQjFxXKPExKMANTkIHRazAVlNyBH4s6d6qUIjIpMstgnI8/T1jcdQNyIeECaRlpWfcLgp+bH0EeQbT1b0vouVDiU28MwFhRNjyqjAzSLl0i5dw5Us+FkRJmvKaeO0dKWBjpkZFFqif9yhUra6qp7Lg62dOlkR+B+ezBcXawY8eZqyzeeyFL5mAn1PZ1o4GfRw7DU9/PA3fnIv6ci8Btz8CGLyA+V8Tu1ETDO9ZCszDawFRA8gr8mH0HbkxyTIEGJDEt55e5qktVAj0CaerblLtuuYtAj0CCPIKo6VGTmh41cbJ3yrdtKLsNj5lrIYXF5UqPTyD1vGk0Mo3HuXOknjtH6vnzqNTU6zfb2+MYEIBjrSA8u3fHsVYtnGoF4VirNk61gjg18H7SLly4oY3KHtFYU3YUtAdnQOtAEpLTOHU5gROX4zgREZ91rD0aQVrG9dFHTW8X6mczOplGKN/9OfGX85Zb0Au2UhsYEekNTATsgWlKqU8t3Ubo9PdxnDIfv5h0NnnbkzpyMN1GvGORujNUBklpSVxLu0ZiWiKJaYkkpSXx+Y7Pabs/gUdDFb6xEOkFP3VL4D/qP3y09aOs0A+ZeDp6EugZyC1et9AxsGMOAxLoEVisBfi+9frivm4XjlPm4xOTTrS3Pakj+9HNyvtTcq+FpF24QPi//0PC9u04+PrmGI2kR0XlKGvn6YlTrVo4N26M59134RiUaURq4VijRoH57au/MlpHNNZYlcJy4Lg7O9AiyJsWQd45yqWkZfBPVEIOo3Picjxzt/9DUur1tZ0qbo45RjqZ51Vda+CWeOOU9zXXGljKJafSGhgRsQf+B9wNhAE7RGSJUuqwpdoInf4+PhPm4mw+EFeNSSd5wlzmX7tMzfsfyjIKianma3qu97mOpLSknO/Tb1xYBuh0KJ1nlitczJksv1h4ZrkCkql1/yCCPIMI9AjMMiDezt551lMQSiljrlYpyMgApVBKEfP779T4ehEqKT2rzzLpV65k1MSjU0dUSgoZKSmolFRUSjIqJSXryEhJQSWn5JAZ8sz7Uq/Lkw1ZRqohTz5+HNLScuqYkkLMLwuyRiFOtWvhctddONYKwqlWrSxDYu9d/P5nkj2isbVyZmg0Rc2Bkx0nBzsaVPekQfWcm5MzMhTnoxM5cTmek9mMz4qDF4m+dn303t9uIJ84TsNNUrJk15QTn6c+xNhS9eY6lXaRX0TuAMYqpXqZ78cAKKU+yev+kizyb7otmKox6TfI0+zgYrZkipLtTywiCIKd2GGHmO+NczsRROzM68Z9mfdePxe4eBn7PJxP0gWcq/ubxiDDWLvLyMhhIHIbjPyuURbfC0dH7BwdEScnxNnZeM06HLFzui6LX7cu7zpEaLJ/X4GjEEuhF7xvDiprn5VSRCakZBmc/yw+yH12G3nTYT41JZILypfP0wazNCOE058WfUbiZl3kDwSy74IMA27LfoOIjARGAvj7+xMaGlqsBvzyMC4A9hngWrsxgh322BvGROwMM5LDnTWna6vK09NVyGFLRHC5kPfcqZ2CuPr1USLGQl7WAYgdCMY1BOzMa5nniFku8948yongsfg38lJTATHPPoNycEA5OICDI8ox5znmNeXoCPb2xYoyXW3fPuxzTX0BpFepwvpNm4pcT2mIj48v9nekoqP7XDkJAnxdhCVJISxJCclxzddFLNb/ymxg8vsdvP5GqSnAFDBGMMV9atnkbZ/nCOaqtz0hcxcXq67icLx7jzwXnh1r1qTZjOlWaxfg+PYd+bZ9x+jRVms3Zsxbea6F1BrzFsFl9LRZWZ9sC0L3ufLytvf5PJ0L3u7fgm4W2uBZ8XfC5U8YkD10bRBw4y9jKUgdOZjkXDMzyY6G3JpUf2U04pIz1ERZLTzbqm3vfv0I+OB9HGrWBBEcatYk4IP39VqIRlNCBrQO5JP7WxDo4woYGzstHT2gMo9gdgANRaQucB54GHjUkg10G/EOoZDLo8pyXmT5YctUurZc9NbZHTUay1LUAJ8lpdIaGKVUmoi8CKzCcFOeoZQ6ZOl2uo14B0a8U+bDalv+2Ga2fbNMJWg0mpJRaQ0MgFJqObDc1npoNBrNzUhlXoPRaDQajQ3RBkaj0Wg0VkEbGI1Go9FYBW1gNBqNRmMVKm2omOIiIpeBs6Woohpws8Vwv9n6fLP1F3SfbxZK0+dblFJ+eV3QBsZCiMjO/OLxVFZutj7fbP0F3eebBWv1WU+RaTQajcYqaAOj0Wg0GqugDYzlmGJrBWzAzdbnm62/oPt8s2CVPus1GI1Go9FYBT2C0Wg0Go1V0AZGo9FoNFZBG5hSIiK9ReSYiJwQkbdsrY+1EZFaIrJORI6IyCERGWVrncoKEbEXkT0i8rutdSkLRMRHRBaIyFHz877D1jpZGxF5xfxeHxSRuSLiUnipioWIzBCRCBE5mE1WVUTWiMhx87VKQXUUFW1gSoGI2AP/A/oAzYBHRKSZbbWyOmnAa0qppsDtwAs3QZ8zGQUcsbUSZchEYKVSqglwK5W87yISCLwMtFNKBWOk+XjYtlpZhVlA71yyt4A/lVINgT/N96VGG5jS0QE4oZQ6pZRKAeYB/W2sk1VRSoUrpXab53EYPzqWS4FXThGRIKAvMM3WupQFIuIFdAGmAyilUpRS0TZVqmxwAFxFxAFww8JZcMsDSqm/gKhc4v7AbPN8NjDAEm1pA1M6AoFz2d6HcRP82GYiInWA1sA2G6tSFkwA3gQybKxHWVEPuAzMNKcFp4mIu62VsiZKqfPAF8A/QDgQo5RabVutygx/pVQ4GA+RQHVLVKoNTOmQPGQ3hd+3iHgAC4HRSqlYW+tjTUTkXiBCKbXL1rqUIQ5AG2CyUqo1kICFpk3KK+a6Q3+gLlATcBeRx22rVcVGG5jSEQbUyvY+iEo4pM6NiDhiGJcflVK/2lqfMqATcJ+InMGYBu0uIj/YViWrEwaEKaUyR6cLMAxOZeYu4LRS6rJSKhX4FehoY53KiksiEgBgvkZYolJtYErHDqChiNQVESeMBcElNtbJqoiIYMzLH1FKfWVrfcoCpdQYpVSQUqoOxme8VilVqZ9slVIXgXMi0tgU9QAO21ClsuAf4HYRcTO/5z2o5I4N2VgCDDXPhwK/WaJSB0tUcrOilEoTkReBVRgeJzOUUodsrJa16QQ8ARwQkb2m7P+UUsttp5LGSrwE/Gg+PJ0CnrSxPlZFKbVNRBYAuzG8JfdQCcPGiMhcoBtQTUTCgHeBT4H5IjICw9A+aJG2dKgYjUaj0VgDPUWm0Wg0GqugDYxGo9ForII2MBqNRqOxCtrAaDQajcYqaAOj0Wg0GqugDYxGUwgiEioi7Qq5p52ITCrkHh8Red6y2hWN3G2LSE3TJVejsRrawGg0FkAptVMp9XIht/kAxTIwYmCJ/9McbSulLiilBlmgXo0mX7SB0WgwAneaOU+mmvlAVouIa7ZbHheRzWaekA55lO+WmSdGRMaaOTdCReSUiGQank+B+iKyV0TGmfe+ISI7RGS/iLyXS5dvMDb91RKRWWbbB0TkFfO++iKyUkR2icgGEWliyv1FZJGI7DOPjrnbNts4aN7vIiIzzbr3iMidpnyYiPxqtnFcRD435fZ56aPR5Ebv5NdortMQeEQp9bSIzAceADJjjrkrpTqKSBdgBhBcSF1NgDsBT+CYiEzGCBYZrJRqBSAiPc02O2AETl1i1v8P0Bh4Uin1vIi0BQLNHCWIiI/ZxhTgWaXUcRG5DfgG6A5MAtYrpQaaOYs88mi7TjZdXwBQSrUwjdRqEWlkXmuFETE72ezH1xiRdvPSR6PJgTYwGs11Tiul9prnu4A62a7NBSOXhoh4iYhPIflRlimlkoFkEYkA/PO4p6d57DHfe2AYnH+As0qprab8FFDP/HFfhmEAPDACMf5ihM0CwNl87Q4MMfVNB2Kk4AyFIcDX5v1HReQskGlg/lRKxQCIyGHgFuBQbn0KqFtzE6MNjEZzneRs5+lA9imy3DGVCouxlLuuvP7XBPhEKfVdDqExukjIakipqyJyK9ALY7QxGBgNRGeOSEpJXmknMrmhH/noM9wCemgqGXoNRqMpGg8BiEgIRiKqmBLUEYcxZZbJKmC4ORpBRAJF5IZETyJSDbBTSi0E3gbamDl4TovIg+Y9Yv7og5Hy9jlTbi9GdsrcbWfnL+Ax8/5GQG3gWH6dyEufonRec/OhRzAaTdG4KiKbAS9K+LSulIoUkU3m4voKpdQbItIU2GJOc8UDj2OMFLITiJFZMvOBcIz5+hgwWUT+Azhi5KrZB4wCppiRcdOB55RSW7K3DfwvW/3fAN+KyAGMKMLDlFLJ2abecpOfPhpNDnQ0ZY1Go9FYBT1FptFoNBqroA2MRqPRaKyCNjAajUajsQrawGg0Go3GKmgDo9FoNBqroA2MRqPRaKyCNjAajUajsQr/D8G6BJvTNYg/AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# collect per-query intersections\n",
    "\n",
    "all_ninters = {}\n",
    "for nprobe in 1, 4, 16, 64: \n",
    "    I = res_per_nprobe[nprobe]    \n",
    "    ninters = [\n",
    "        faiss.eval_intersection(I[i0 : i0 + 1], Igt[i0 : i0 + 1])\n",
    "        for i0 in range(10**6)\n",
    "    ]\n",
    "    all_ninters[nprobe] = ninters\n",
    "    pyplot.plot(np.bincount(ninters), 'o-', label=f\"nprobe {nprobe}\")\n",
    "pyplot.xlabel(\"nb intersections\")\n",
    "pyplot.ylabel(\"nb queries\")\n",
    "pyplot.legend()\n",
    "pyplot.grid()    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "id": "unavailable-forward",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "nprobe 1\n",
      "n=   1000 stddev 0.705 %\n",
      "n=  10000 stddev 0.223 %\n",
      "n= 100000 stddev 0.071 %\n",
      "n=1000000 stddev 0.022 %\n",
      "nprobe 4\n",
      "n=   1000 stddev 0.727 %\n",
      "n=  10000 stddev 0.230 %\n",
      "n= 100000 stddev 0.073 %\n",
      "n=1000000 stddev 0.023 %\n",
      "nprobe 16\n",
      "n=   1000 stddev 0.525 %\n",
      "n=  10000 stddev 0.166 %\n",
      "n= 100000 stddev 0.053 %\n",
      "n=1000000 stddev 0.017 %\n",
      "nprobe 64\n",
      "n=   1000 stddev 0.382 %\n",
      "n=  10000 stddev 0.121 %\n",
      "n= 100000 stddev 0.038 %\n",
      "n=1000000 stddev 0.012 %\n"
     ]
    }
   ],
   "source": [
    "# do some math to compute standard deviations\n",
    "\n",
    "for nprobe in 1, 4, 16, 64: \n",
    "    intersection_measures = np.array(all_ninters[nprobe]) / 10 \n",
    "    variance = intersection_measures.var()\n",
    "    print(\"nprobe\", nprobe)\n",
    "    for n in 10**3, 10**4, 10**5, 10**6: \n",
    "        # sum of independent variables\n",
    "        # https://en.wikipedia.org/wiki/Variance#Sum_of_uncorrelated_variables_.28Bienaym.C3.A9_formula.29\n",
    "        variance_of_sum = n * variance\n",
    "        variance_of_mean = variance_of_sum * (1 / n) ** 2\n",
    "        sttdev_of_mean = np.sqrt(variance_of_mean)\n",
    "        print(f\"n={n:-7} stddev {100*sttdev_of_mean:.3f} % (percentage points)\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "realistic-guinea",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
